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ABSTRACT 

FACTUNC  is  a system  for  solving  unconstrained  minimization  problems  based 
on  the  concept  of  factorable  programming.  This  concept  enables  the  user  to 
provide  the  problem  function  and  data  in  a user  friendly  way  and  does  not 
require  user-supplied  derivatives.  The  system  utilizes  the  factorable 
function  concept  to  obtain  the  first  and  second  derivatives  required  for 
unconstrained  optimization.  In  all  cases  derivatives  are  obtained  rapidly  and 
accurately  (up  to  roundoff  errors  due  to  machine  precision) , as  opposed  to 
finite  differencing. 

As  a system  for  nonlinear  minimization,  FACTUNC  allows  several  options. 
First  the  user  can  solve  regression  (nonlinear  least  squares)  problems  by 
providing  the  regression  equation  and  the  data  for  the  dependent  and 
independent  variables.  The  second  option  allows  for  the  minimization  of  the 
sum  of  an  indexed  function.  The  user  provides  the  function,  and  the  indexed 
data.  This  can  be  used  for  example  to  solve  maximum  likelihood  estimation 
problems  when  the  user  provides  the  negative  of  the  (weighted)  log  of  the 
frequency  function  and  the  data.  The  third  option  is  simply  to  minimize  a 
function  supplied  by  the  user.  Utilizing  barrier  function  methodology,  this 
third  option  can  sometimes  be  used  to  solve  constrained  problems. 
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1 . INTRODUCTION 


We  report  here  on  FACTUNC,  a system  for  solving  unconstrained  minimization 
problems  developed  at  the  National  Bureau  of  Standards.  FACTUNC  differs  from 
other  optimization  systems  in  one  fundamental  aspect:  the  internal 
representation  of  functions  as  "factorable”.  In  section  2 we  discuss  the 
concept  and  merits  of  factorable  functions  in  more  detail.  One  advantage  of 
this  approach  is  that  once  a function  is  represented  in  factorable  form,  its 
gradient,  Hessian  and  even  higher  order  derivatives  may  be  calculated 
automatically  and  exactly.  In  order  to  use  FACTUNC,  the  user  need  only  supply 
the  system  with  a statement  of  the  problem  function  in  a FORTRAN- like 
expression.  The  system  then  proceeds  to  translate  this  function  statement 
into  factorable  form,  and  uses  this  representation  to  compute  the  gradient  and 
Hessians  as  required. 

The  Hessian  matrices  computed  by  FACTUNC,  are  initially  obtained  in  dyadic 
form,  i.e.,  as  a sum  of  rank  one  matrices.  Moreover,  the  factors  of  these 
rank  one  matrices  are  obtained  from  the  gradient  evaluation.  The  advantages 
of  this  representation,  and  how  it  can  be  utilized  are  explained  in  more 
detail  in  section  2. 

There  are  a number  of  different  ways  to  use  FACTUNC.  For  example,  FACTUNC 
may  be  used  as  a general  unconstrained  minimization  system.  The  system  also 
provides  a user-friendly  front-end  for  objective  functions  which  involve  an 
indexed  sum  of  terms  of  the  same  form  (such  as  the  negative  of  the  logarithm 
of  the  likelihood  function  which  is  used  to  solve  maximum  likelihood) . When 
working  with  this  option,  the  user  need  only  supply  the  general  form  of  a 
typical  term  in  the  indexed  sum  (e.g.,  a typical  term  in  the  logarithm  of  the 
likelihood  function) . The  nonlinear  least  squares  problem  is  another  example 
of  an  indexed  sum.  Since  nonlinear  regression  problems  arise  frequently  in 
applications,  FACTUNC  has  a special  option  for  solving  this  type  of  problem. 

In  order  to  solve  such  a least  squares  minimization  problem,  the  user  need 
only  supply  the  general  form  of  the  nonlinear  regression  equation,  and,  of 
course,  the  data.  Finally,  by  using  a sequential  unconstrained  minimization 
approach  for  solving  constrained  minimization,  FACTUNC  may  be  incorporated  as 
part  of  a constrained  optimization  system.  Thus,  FACTUNC  may  be  used  to  solve 
the  sequence  of  unconstrained  problems  which  arise  in  the  barrier  method,  or 
in  the  method  of  centers.  In  each  of  these  cases,  the  option  of  using  an 
indexed  sum  can  be  invoked. 

The  use  of  the  various  options  available  in  FACTUNC  is  explained  in  detail 
in  sections  3-5.  For  ease  of  presentation,  we  have  chosen  to  start  with  the 
nonlinear  least  squares  option.  The  underlying  approach  in  implementing 
FACTUNC  is  common  to  all  the  various  options,  and  is  described  in  section  1.1 
below. 
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1 . 1 Representation  of  Functions  in  FACTUNC 

Computer  programs  for  solving  nonlinear  problems  are  widely  available,  A 
major  difficulty  in  using  the  software  is  addressed  in  the  FACTUNC  system. 

This  difficulty  is  how  to  construct  a user-friendly  way  to  represent  the 
information  which  is  required  for  the  nonlinear  minimization,  concerning  the 
nonlinear  function,  which  can  also  be  accessed  in  a computationally  efficient 
manner  by  algorithms  for  minimization.  It  is  solved  in  FACTUNC  by  allowing 
the  user  to  provide  a FORTRAN- like  expression  to  represent  the  nonlinear 
objective  function.  This  is  then  processed  by  the  system  and  represented 
internally  in  the  form  of  a factorable  function.  In  this  form  the  quantities 
required  by  algorithms  for  solving  the  problem:  function  values,  first 
derivatives,  and  second  derivatives  are  automatically  computed  by  the  system. 

The  statement  for  the  function  allows  the  regular  operations  such  as  plus 
(+) , minus  (-),  product  (*) , division  (/) , power  (**) , and  parentheses.  If 
parentheses  do  not  conform  (e.g.,  an  extra  left  parenthesis,  or  some  other 
mismatch,  an  error  statement  is  produced).  In  addition,  the  statement  allows 
for  single  variable  functions.  The  single-variable  functions  available  for 
use  by  the  system  are  given  in  Table  1.  Note  that  any  single-variable 
function  can  theoretically  be  used.  These  are  just  the  ones  programmed  to 
date . 

FACTUNC  also  accepts  complex  terms  which  involve  the  nesting  of  parentheses 
and  single  variable  transformations  within  other,  provided  that  the  resulting 
expression  represents  a well  defined  mathematical  expression.  Thus  FACTUNC 
will  be  able  to  accept  and  process  a function  of  the  form 

A*SQRT(X**3+  (Y/(X+Y))*ENTR0P(Z))-1.7*DECAY(SIN(Y)) . 

However,  the  expression 


(GAMMA(X+Y) )/(2*L0G(Y) ) ) 

is  incorrect.  The  system  will  notify  the  user  of  the  excess  in  the  number  of 
right  parentheses.  Removal  of  the  last  parenthesis,  will  result  in  an 
expression  acceptable  to  FACTUNC. 


TABLE  1.  SINGLE  ARGUMENT  FUNCTIONS  AVAILABLE  IN  THE  FACTUNC  SYSTEM 


FACTUNC  NAME 

ALGEBRAIC  NOTATION 

DESCRIPTION 

EXP(X) 

e'* 

Exponential 

LOG(X) 

logioX 

Common  Logarithm 

LN(X) 

log^x 

Natural  Logarithm 

SIN(X) 

sin  X 

Sine 

COS(X) 

cos  X 

Cosine 

TAN(X) 

tan  X 

Tangent 

GAMMA(X) 

r(x) 

Gamma 

ARCSIN(X) 

arcs in  x 

Arcsine 

ARCTAN(X) 

arc tan  x 

Arctangent 

CUMNOR(X) 

1 

— EXP(-t2/2)dt 

^ J« 

Cumulative  Normal 

SQRT(X) 

Jx 

Square  Root 

DECAY(X) 

1-e'* 

Decay  Function 

ENTROP(X) 

X'in  X 

Entropy 

NORDEN(X) 

— EXP(-x2/2) 

J 2ir 

Normal  Density 

BARLN(X) 

< 

® , x<0 

-In  X,  x>0 

Logarithmic 

Barrier  Function 
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2.  MATHEMATICAL  BACKGROUND  ON  FACTORABLE  FUNCTIONS 

This  section  is  provided  as  background  for  the  interested  reader.  It  is 
not  necessary  to  understand  the  theory  of  factorable  functions  in  order  to  use 
the  FACTUNC  system.  Those  interested  only  in  using  the  system  may  therefore 
skip  this  section. 

A factorable  function  is  a multivariable  function  that  can  be  written  as 
the  last  of  a finite  sequence  of  functions,  in  which  the  first  n functions  in 
the  sequence  are  just  the  coordinate  variables,  and  each  function  beyond  the 
nth  is  a sum,  a product,  or  a single-variable  transformation  of  previous 
functions  in  the  sequence.  More  rigorously,  let  [f^(x),  f2 (x) , ...,  fL (x) ] be 
a finite  sequence  of  functions  such  that  f^  :R"  -*■  R,  where  each  f^  (x)  is 
defined  according  to  one  of  the  following  rules. 

Rule  1 ■ For  i=l , ...,  n,  f^ (x)  is  the  value  of  the  i^^  Euclidean  coordinate: 

fi(x)=Xi.  (2.1) 

Rule  2 . For  i=n+l , ...,  L,  f^ (x)  is  formed  using  one  of  the  following 

compositions : 


a. ) 

fi  (x) 

= or 

b.) 

fi  (x) 

= fj(i,(x)  * fk(i,(x);  or 

(2.2) 

c. ) 

fi  (x) 

~ ( l ) 1 > 

where  j(i)  < i,  k(i)  < i. 

and  Tl 

is  a function  of  a single  variable. 

Then 

f(x)  = fL  (x)  is  a factorable  function  and  [fj^(x),  f2  (x) , ...,  fL  (x)  ] is  a 
factored  sequence.  Thus  a function,  f(x),  will  be  called  factorable  if  it  can 
be  formed  according  to  Rules  1 and  2,  and  the  resulting  sequence  of  functions 
will  be  called  a factored  sequence  or,  at  times,  the  function  written  in 
factored  form. 

Although  it  is  not  always  immediately  grasped,  the  concept  of  a factorable 
function  is  actually  a very  natural  one.  In  fact,  it  is  just  a formalization 
of  the  natural  procedure  one  follows  in  evaluating  a complicated  function. 
Consider  for  example  the  fuhction 

f(x)  = (a^x) (sin[b^x] ) (exp [c^x] ) , (2.3) 

where  a,b,c  and  x are  (2x1)  vectors.  The  natural  approach  to  evaluating 
this  function  for  specific  values  of  and  X2  is  first  to  compute  the 
quantities  within  the  parentheses,  then  to  apply  the  sine  and  exponential 
functions,  and  finally  to  multiply  the  three  resulting  quantities.  This  might 
be  done  in  stages  as  follows. 


6 


(2.4) 


fe 

= b,f, 

fii 

■“  f 9 f 1 0 

f7 

- b2f2 

fl2 

= sin(f8) 

= a^fi 

fs 

= fg  4-  f? 

fia 

= exp ( f 1 1 ) 

” ^2^2 

fg 

= c^fi 

fl  4 

- fs  • fl2 

- fa  + f4 

fio 

= C2f2 

fl5 

" fia  * fi4 

This  is  one  possible  factored  sequence  for  the  function  in  (2.3). 


In  order  to  appreciate  fully  the  value  of  factorable  functions  the  concept 
of  an  outer  product  matrix  must  be  introduced.  An  (m  x n)  matrix  A is  called 
an  outer  product  matrix  if  there  exists  a scalar  a,  an  (m  x 1)  vector  a,  and 
an  (n  X 1)  vector  b such  that 


A = aab^  . 

The  expression  aab^  is  called  an  outer  product  or  a dvad.  Note  that  a dyad  is 
conformable  since  the  dimensions  of  the  product  are  (m  x 1)(1  x 1)(1  x n) , 
which  yields  the  (m  x n)  outer  product  matrix  A as  desired.  A useful  property 
of  outer  product  matrices  is  that,  if  they  are  kept  as  dyads,  matrix 
multiplication  with  them  is  simplified  to  inner  products  alone,  saving  the 
computations  required  to  form  the  matrices  involved.  For  example, 

Ac  = aQ:[b^c]  , 

d"^  A - [ d^  a ] ab^  , and 

AF  = aa[b^F] , 

where  c is  (n  x 1),  d is  (m  x 1)  and  F is  (n  x m) . 

Factorable  functions  possess  two  very  special  properties  that  can  be 
exploited  to  produce  efficient  (fast  and  accurate)  algorithms:  (i)  once 
written  in  factorable  form,  their  gradients  and  Hessians  may  be  computed 
exactly,  automatically,  and  efficiently;  and  (ii)  their  Hessians  occur 
naturally  as  sums  of  dyads  whose  vector  factors  are  gradients  of  terms  in  the 
factored  sequence.  The  first  of  these  properties  is  utilized  in  FACTUNC  for 
providing  the  first  and  second  derivatives  of  nonlinear  functions.  The  second 
has  obviated  the  task  of  multiplying  a matrix  by  a vector,  reducing  it  to  a 
series  of  inner  products,  as  noted  above,  which  in  many  cases  results  in  less 
effort. 

In  order  to  clarify  these  concepts,  we  consider  again  the  illustrative 
function  in  (2.3).  Table  1 is  a display  of  the  gradient  and  Hessian  of  this 
function.  The  entries  in  each  column  are  the  summands  in  the  expressions  for 
the  gradient  and  Hessian.  For  example. 
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Vf  = [sin(b^x) ] [exp(c^x) ] a + [a^x] [cos(b^x) ] [exp(c^x) ]b 
+ [a’-x]  [sin(b^x)  ] [exp(c'^x)  ]c. 

The  table  also  illustrates  the  left- to-right , tree-like  structure  of  the 
derivatives  involved.  From  the  table  it  can  be  seen  that  both  the  gradient 
and  Hessian  naturally  have  the  dyadic  structure  discussed  above.  Notice  too 
that  the  vectors  in  the  monads  and  dyads  are  drawn  from  the  set  {a,b,c),  each 
of  which  is  the  gradient  of  a factored  sequence  function  in  (2.4). 


TABLE  1. 

MONADIC  AND  DYADIC  STRUCTURE  OF  GRADIENT 
AND  HESSIAN  OF  ILLUSTRATIVE  FUNCTION 


Function 

Gradient  Summands 

Hessian  Summands 

a^  X • s in  [ b"^  X ] exp  [ c*^  x ] 

(sin[b^x] exp [c^x] :a) 

(cos [b^x] exp [c^x] :ab) 

(sin[b’-x] exp [c^x]  :ac) 

(a^x'cos [b^x] exp [c^x] :b) 

(cos [b^x] exp [c^x] :ba) 

(-a^x*sin[b^x]exp[c^x] :bb) 

(a^x«cos [b^x] exp [c'^x]  :bc) 

(a^x*sin[b^x] exp [c^x] :c) 

(sin[b^x] exp[c^x] :ca) 

(a'^x’cos  [b^x] exp [c'^x]  :cb) 

(a^x« sin[b^x] exp [c^x] :cc) 

It  is  important  to  understand  that  the  derivative  calculations  performed 
by  the  FACTUNC  system  are  not  estimations,  but  mathematically  exact 
calculations.  Furthermore,  they  are  also  compact,  since  factored  sequences 
mimic  hand  calculations , and  thus  this  technique  is  different  from  symbolic 
manipulation  techniques  for  differentiation,  which  tend  to  produce  large 
amounts  of  code.  The  techniques  used  in  Factorable  Programming  are  efficient 
exploitations  of  the  special  structure  inherent  in  factorable  functions  and 
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their  partial  derivative  arrays.  These  results  were  extended  to  higher  order 
derivatives  by  Jackson  and  McCormick  [1986],  who  showed  that  factorable 
functions  have  arrays  of  n'^^  order  derivatives  (tensors)  which  are  naturally 
computed  as  sums  of  generalized  outer  product  matrices  (polyads) . 


3.  SOLVING  REGRESSION  (LEAST  SQUARES)  PROBLEMS  (OPTION  1) 

We  begin  the  discussion  of  the  use  of  the  code  with  a description  of  the 
simplest  way  to  use  it;  for  solving  nonlinear  regression  problems. 

A regression  equation  relates  a dependent  variable  to  a function  involving 
one  or  more  independent  variables  and  parameters.  Given  a set  of  values  of 
the  dependent  variable  and  simultaneously  observed  values  of  the  independent 
variables,  the  "best"  values  of  the  parameters  in  a least  squares  sense  are 
those  which  minimize  the  sum  of  the  squared  differences  between  the  observed 
dependent  values  and  the  estimated  ones.  Models  giving  rise  to  least  squares 
problems  are  legion  in  the  physical  and  social  sciences.  The  FACTUNC  system 
provides  a user  friendly  way  to  represent  the  data  and  the  regression 
function.  This  is  particularly  important  when  the  parameters  enter 
nonlinearly  in  the  regression  function. 

In  the  FACTUNC  system,  one  variable,  y,  is  the  dependent  variable  and 
there  may  be  one  or  more  independent  variables  , . . . . The  problem  modeler 

must  specify  the  form  of  the  function  f which  is  believed  to  represent  the 
relationship  between  the  variables.  The  function  has  parameters  a^^  , . . . , ajjj 
which  must  be  calculated  in  a way  to  give  a "best"  fit  to  the  observed  data. 
Ideally,  the  relationship  y = f(aj^  , . . . ,an,;xj^  , . . . ,x„)  would  hold  precisely, 
i.e.  , 

y^  = f(ai , . . . , . . . , x^) 

for  each  i=l,2,...,k,  where  k is  the  total  number  of  observations  and  the 
super/subscript  i indicates  the  values  of  the  variables  at  the  i^^ 
observation.  However,  it  will  not  usually  be  possible  to  get  an  exact  fit. 

In  the  least  squares  method,  the  parameters  are  chosen  to  minimize  the 
quantity 

k 

S [y^  -f(ai  , . . . , . . . , x^)]^  , 

i=l 

the  sum  of  the  squares  of  the  residuals.  A least  squares  problem  is  said  to 
be  nonlinear  if  the  parameters  appear  nonlinearly  in  the  regression  equation. 
For  more  discussion  see  McCormick  [1983],  pp.  93-99. 

To  solve  a nonlinear  least  squares  problem  using  FACTUNC,  the  user  must 
provide  an  input  file  giving  the  regression  equation,  the  names  of  the 
parameters  along  with  their  starting  values  and  upper  and  lower  bounds,  and  a 
list  of  the  observed  data.  The  general  structure  of  the  input  file  is 
explained  in  the  following. 
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Input  File  Format 


( 1 ) NLSQ 

All  characters  on  this  line  are  ignored  and  can  be  used  as  a comment 
to  entitle  the  problem. 

(2)  This  can  also  be  used  to  identify  the  regression  equation. 

(3)  Regression  equation  - given  in  the  following  format: 

(dependent  variable)  = (expression  involving  functions, 
parameters,  constants,  and  independent  variabies)$ 

The  notation  follows  the  usual  FORTRAN  conventions  (e.g. 
indicates  multiplication) . See  the  table  for  a list  of  the  allowed 
functions.  The  equation  may  be  written  on  several  lines,  but  the  ”=" 
must  appear  on  the  first  line.  The  end  of  the  equation  is  signalled 
by  the  "$".  Spaces  are  ignored. 

(4)  Comment  - demarcates  problem  constants  section  (must  be  present) . 

(5)  Problem  constants  - number  of  parameters,  constants,  independent 
variables,  and  observations  (in  free  format). 

(6)  Comment  - used  to  demarcate  parameter  section  (must  be  present) . 

(7)  Parameter  section  - A line  is  given  for  each  parameter  in  the 
regression  equation  showing  the  parameter  name  in  quotes,  a starting 
value,  a lower  bound,  and  an  upper  bound  (in  free  format). 

(8)  Comment  - used  to  demarcate  constant  section  (must  be  present) . 

(9)  Constant  section  - A line  is  given  for  each  named  constant  in  the 
regression  equation  showing  the  constant  name  in  quotes  and  its 
value.  This  section  may  be  blank,  but  the  preceding  comment  (8)  must 
always  be  present. 

(10)  Comment  - used  to  demarcate  data  section  (must  be  present). 

(11)  Data  section  - A line  is  given  showing  the  names  of  the  variables  (in 
quotes)  in  the  regression  equation.  Under  each  variable  name  a 
column  of  data  values  is  listed.  Each  row  of  this  section  represents 
one  observation. 


As  an  example  of  the  use  of  the  regression  (least  squares)  option  consider 
the  data  in  Table  2.  The  elements  in  the  table  are  the  ratios  of  those  people 
who  died  in  the  category  indicated  in  the  appropriate  row  and  column,  divided 
by  those  who  would  be  expected  to  die  if  the  sample  were  taken  from  a general 
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population  of  nonsmokers.  If  the  matrix  elements  are  denoted  by  DEATH,  the 
row  indicator  by  DEPTH,  and  the  column  indicator  by  GIGS,  a reasonable 
regression  equation  is 

DEATH  = 1 + (a^  + a2 DEPTH) (1  - EXP( -ag GIGS) ) 

where  (a^  , a2  , are  parameters.  The  values  of  DEPTH  and  GIGS  need  to  be 

quantified.  The  values  used  for  this  example  are  given  in  Figure  1,  the  input 
file  for  FAGTUNG.  The  optimization  problem  is: 

16 

min  S [DEATHi  - {1  -f  (a^  + a2DEPTHi)(l  - EXP( -ag  GIGS^  ) ) ) ] ^ 
a^^  i“l 


TABLE  2.  DATA  ON  SMOKING  AND  HEALTH 


Degree 

of 

Inhalation 

Number  of  Gigarettes  per  Day 

1-9 

10-19 

20-40 

40+ 

None 

1.29 

1.46 

1.56 

2.05 

Slight 

1.29 

1.68 

1.84 

1.97 

Moderate 

1.61 

1.82 

1.84 

2.01 

Deep 

1.88 

1.76 

2.18 

2.50 

The  general  format  of  the  input  file  for  this  example  is  given  in  Figure  1, 
with  each  line  of  the  file  numbered  and  keyed  to  the  input  file  format  given 
above.  Figure  2 is  the  exact  input  file  for  this  example,  while  Figure  3 is 
part  of  the  FAGTUNG  output.  Table  1 is  a list  of  the  functions  of  a single 
argument  which  can  be  used  in  describing  input  functions. 
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FIGURE  1.  SAMPLE  INPUT  FILE  FOR  LEAST  SQUARES  OPTIONS. 


(1) 

NLSQ  - EXCESS  DEATHS 

AS  FCN  OF  DAILY  SMOKES  AND  DEPTH  OF  INHALATION 

(2) 

REGRESSION 

EQUATION 

(3) 

DEATH=1 . + ( A1*DEPTH+A2 ) * ( C 1 - EXP ( - A3*CIGS ) ) $ 

(4) 

NUMBER  OF  PARAMETERS, 

CONSTANTS,  INDEPENDENT  VARIABLES,  OBSERVATIONS 

(5.) 

3 1 

2 16 

(6) 

PARAMETERS , 

STARTING 

VALUE,  LOWER  BOUND,  UPPER  BOUND 

(7) 

'Al' 

0.1 

0.  10. 

'A2' 

0.6 

0.  10. 

'A3' 

0.01 

0.  10. 

(8) 

CONSTANTS : 

NAME,  VALUE 

(9) 

'Cl' 

1. 

(10) 

VARIABLES : 

NAME  ABOVE  COLUMN 

(11) 

'DEATH' 

'CIGS' 

'DEPTH' 

1.29 

5. 

0. 

1.46 

15. 

0. 

1.56 

30. 

0. 

2.05 

65. 

0. 

1.29 

5. 

.3333 

1.68 

15. 

.3333 

1.84 

30. 

.3333 

1.97 

65. 

.3333 

1.61 

5. 

.6667 

1.82 

15. 

.6667 

1.84 

30. 

.6667 

2.01 

65. 

.6667 

1.88 

5. 

1.0000 

1.76 

15. 

1.0000 

2.18 

30. 

1.0000 

2.50 

65. 

1.0000 

*A11  inputs  are  free  form.  The  numbers  in  parentheses  correspond  to  the 
numbers  in  the  explanation  above,  but  are  not  part  of  the  input  file. 
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FIGURE  2.  Exact  Input  File  for  Smoking  and  Health  Regression  Problem 


NLSQ  - SMOKING  AND  HEALTH  REGRESSION 
REGRESSION  EQUATIONS 

DEATH  = 1.  +(A1*DEPTH+A2)*(C1-EXP(-A3*CIGS))$ 

NUMBER  OF  PARAMETERS , CONSTANTS , INDEPENDENT  VARIABLES , OBSERVATIONS 
3 2 2 16 

PARAMETERS  , STARTING  VALUE,  LOWER  BOUND,  UPPER  BOUND 
'Al'  4.  0.  10. 

'A2'  3.  0.  10. 

'A3'  .1  0.  10. 

CONSTANTS:  NAME,  VALUE 
'DUM'  2.2 

'Cl'  1. 

VARIABLES:  NAME  ABOVE  COLUMN:  THOSE  NOT  HERE  ARE  ASSUMED  ON  FILES 
'DEATH'  'GIGS'  'DEPTH' 


1.29 

1.46 

1.56 

2.05 

1.29 

1.68 

1.84 

1.97 

1.61 

1.82 

1.84 

2.01 

1.88 

1.76 

2.18 

2.50 


5. 

15. 

30. 

65. 

5. 

15. 

30. 

65. 

5. 

15. 

30. 

65. 

5. 

15. 

30. 

65. 


0. 

0. 

0. 

0. 

3333 

3333 

3333 

3333 

6667 

6667 

6667 

6667 


1.000 

1.000 

1.000 

1.000 
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FIGURE  3.  Partial  Output  for  Least  Squares  Exanple. 


NLSQ  - SMOKING  AND  HEALTH  REGRESSION 
LEAST  SQUARES  PROBLEM 
HESSIAN  FORMED  EXPLICITLY 
REGRESSION  EQUATION 

DEATH  = 1.  +(A1*DEPTH+A2)*(C1-EXP(-A3*CIGS))$ 

SINGLE  CHARACTER  OUTPUT  STRING  OF  (FUNCTION)  EQUATION 
1.+(A1*DEPTH+A2)*(C1-EXP(-A3*CIGS))$ 

NUMBER  OF  PARAMETERS , CONSTANTS , INDEPENDENT  VARIABLES , OBSERVATIONS 
THE  NUMBER  OF  UNKNOWNS  IS  = 3 

THE  NUMBER  OF  SYMBOLIC  CONSTANTS  IN  THE  EQUATION  (FUNCTION)  = 2 
THE  NUMBER  OF  INDEPENDENT  VARIABLES  = 2 

THE  NUMBER  OF  DATA  OBSERVATIONS  = 16 


PARAMETERS  , 

STARTING  VALUE, 

LOWER  BOUND,  UPPER 

. BOUND 

Al 

.400000E+01 

.OOOOOOE+00 

lOOOOOE+02 

A2 

.300000E+01 

.OOOOOOE+00 

lOOOOOE+02 

A3 

.lOOOOOE+00 

.OOOOOOE+00 

lOOOOOE+02 

CONSTANTS:  NAME,  VALUE 

DUM  .220000E+01 

Cl  .lOOOOOE+Ol 

VARIABLES:  NAME  ABOVE  COLUMN 
DEATH  CIGS 

DEPTH 

1.29000 

5.00000 

.00000 

1.46000 

15.00000 

.00000 

1.56000 

30.00000 

.00000 

2.05000 

65.00000 

.00000 

1.29000 

5.00000 

.33330 

1.68000 

15.00000 

.33330 

1.84000 

30.00000 

.33330 

1.97000 

65.00000 

.33330 

1.61000 

5.00000 

.66670 

1.82000 

15.00000 

.66670 

1.84000 

30.00000 

.66670 

2.01000 

65.00000 

.66670 

1.88000 

5.00000 

1.00000 

1.76000 

15.00000 

1.00000 

2.18000 

30.00000 

1.00000 

2.50000 

65.00000 

1.00000 

THE  OPTIMAL  VALUES  OF  THE  UNKNOWNS  ARE 

1 A1  = .54488E+00 

2 A2  = .74750E+00 

3 A3  = .93753E-01 

THE  LEAST  SQUARES  VALUE  IS  = .46916E+00 

THE  GRADIENT  IS  -.14019E-09  -.25256E-09  -.23803E-08 
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4.  MINIMIZING  THE  SUM  OF  AN  INDEXED  FUNCTION  (OPTION  2) 


The  general  structure  of  the  input  file  for  minimizing  the  sum  of  an 
indexed  function  is  explained  in  the  following.  This  option  is  triggered  by 
placing  the  word  MAXL  in  columns  1-4  of  the  first  line  in  the  input  file.  The 
main  difference  between  the  input  format  of  this  option  and  the  least  squares 
option  is  that  on  line  three  the  "="  sign  and  name  of  the  dependent  variable 
do  not  appear.  The  user  merely  writes  on  line  3 (and  for  as  many  following 
lines  as  needed)  the  indexed  function  whose  sum  is  to  be  minimized.  The 
remainder  of  the  input  file  is  exactly  that  described  in  Section  3. 

4.1  MAXIMUM  LIKELIHOOD  ESTIMATION 

Another  common  method  of  finding  estimates  of  parameter  values  that 
explain  observed  data  is  the  method  of  maximum  likelihood  estimation.  Let  b = 
(bj^  , . . . ,bj. ) be  r unknown  parameters  of  a frequency  function  go(yfb)  of  the 

random  variable  y.  Let  y^^ yjj,  be  m observations  of  y.  The  likelihood 

function  associated  with  these  observations  and  frequency  function  is 

L(yi  , • . . .y^  , . . . ,bj.)  = go  (Xi  ,b)go  (y2  ,b)  . . .gg  (y„  ,b)  . (4.1) 

An  important  method  of  estimating  the  values  of  bj^ b^  based  on  the 

observed  . ,yj^  is  to  maximize  the  likelihood  function  (4.1)  . The 

bj^ bj.  that  maximize  (4.1)  are  called  maximum  likelihood  estimators  and 

have  many  desirable  statistical  properties.  The  reader  is  referred  to  Cramer 
[??]  for  a fuller  discussion. 

The  problem  of  maximizing  L(yj^ y^, , b)  is  an  unconstrained 

mathematical  programming  problem.  Since  the  logarithm  of  the  likelihood 
function  achieves  its  maximum  at  the  same  b as  the  likelihood  function  itself, 
the  general  problem  of  likelihood  estimation  is  stated  as  follows. 

Find  values  (bj^  , ...,  b^ ) that 

r 

maximize  S In  goCXiib),  (^.2) 

i=l 
or 

r 

minimize  -S  In  go(yi,b),  (4.3) 

i=l 


The  example  for  this  section  comes  from  the  world  of  biomedicine  and  is 
taken  from  [Bracken  and  McCormick,  1968] . The  mathematical  technique 
illustrated  is  the  minimization  of  the  negative  of  the  log- likelihood 
function. 

It  is  hypothesized  that  the  population  of  systolic  blood  pressures  can  be 
separated  into  three  separate  groups.  The  distribution  of  blood  pressures 
within  each  of  these  groups  can  be  represented  by  a normal  frequency  function 
bet  pj^  , P2  , and  P3  represent  the  proportions  of  the  population  in  each  of  the 
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three  groups.  Let  (Ms.^Ts)  he  the  means  and  standard 

deviations  of  the  normal  frequency  functions  corresponding  to  each  group. 
These  nine  values  correspond  to  the  unknovm  {bj } parameters. 


Then  under  these  assumptions  the  frequency  function  for  the  random 
variable  y,  which  denotes  systolic  blood  pressure,  is  obtained  by  summing  the 
frequency  functions  of  the  individual  groups  times  their  probability  of 
occurrence  to  yield 


1 

J 2n 


k=l 


exp 


(y  - 


where 


Pi  + P2  + Pa  = 1- 


There  are  eight  parameters  in  this  frequency  function  since  one 
proportion,  or  probability,  can  be  eliminated.  Let  pg  = l-p3^-p2.  Let  n^ 
equal  the  frequency  of  occurence  of  the  i^^  observation.  The  problem  then  is 
to  find  values  of  (pj^  , P2  , » A*2  » ^3  » » ^2  » ^a  ) that 


minimize 


1 

Pi 

P2 

exp 

^1 

2af 

+ exp 

^2 

I-P1-P2 

exp 

^3 

The  data  for  this  are  given  in  Table  4.  Using  the  FACTUNC  Systems  yields 


estimates , 

. (See 

Appendix  A) 

Pi  = 

.255 

= 128.8 

P2  = 

.594 

= 158.3 

Pa  = 

.151 

^3  = 222.2 

= 10.3 


a2  = 21.7 


03  = 18.5 


The  input  file  for  this  problem  is  given  in  Figure  4.  Abbreviated  output 
is  given  in  the  Appendix  A. 
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TABLE  4.  SYSTOLIC  BLOOD  PRESSURE  VALUES  WITH  FREQUENCY  OF  OCCURRENCE 


Systolic 

Blood 

Pressure 

Frequency 

of 

Occurrence 

Systolic 

Blood 

Pressure 

Frequency 

of 

Occurrence 

Systolic 

Blood 

Pressure 

Frequency 

of 

Occurrence 

95 

1 

150 

17 

200 

3 

105 

1 

155 

4 

205 

3 

110 

4 

160 

20 

210 

8 

115 

4 

165 

8 

215 

1 

120 

15 

170 

17 

220 

6 

125 

15 

175 

8 

225 

0 

130 

15 

180 

6 

230 

5 

135 

13 

185 

6 

235 

1 

140 

21 

190 

7 

240 

7 

145 

12 

195 

4 

245 

1 

260 

2 
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FIGURE  4.  INPUT  FILE  FOR  SYSTOLIC  BLOOD  PRESSURE  PARAMETERS  ESTIMATION. 


MAXL-FIND  PARAMETERS  DEFINING  FREQUENCY  FUNCTION  MIN  -(LOG  LIKE  FN) 
FREQUENCY  FUNCTION  IS  WEIGHTED  SUM  OF  THREE  NORMALS 
( -FREQ) *LN( (P/S IG)*NORDEN( (Y-MU)/SIG)+( ( (1 . -P-Q) /SIG2)* 

NORDEN( (Y-MU2)/SIG2) )+(Q/SIG3)*NORDEN( (Y-MU3)/SIG3) ) $ 

NUMBER  OF  UNKNOWN  PARAMETERS , CONSTANTS , IND . VARS.  AND  OBSERVATIONS 
8 0 2 31 

PARAMETERS:  STARTING  VALUE,  LOWER  AND  UPPER  BOUND 
'MU'  133.2425  120.  200. 

'SIG'  9.47758  8.  70. 

'P'  .28802  .2  .9 

'MU2'  171.10458  141.  230. 

'SIG2'  33.76752  10.  100. 

'Q'  .2  .1  .9 

'MU3'  200.  100.  300. 

'SIG3'  30.  10.  50. 

CONSTANTS  ( NONE  FOR  THIS  PROBLEM) 

INDEXED  DATA 


'Y' 

'FREQ' 

95. 

1. 

105. 

1. 

110. 

4. 

115. 

4. 

120. 

15. 

125. 

15. 

130. 

15. 

135. 

13. 

140. 

21. 

145. 

CM 

150. 

17. 

155. 

4. 

160. 

20. 

165. 

8. 

170. 

17. 

175. 

8. 

180. 

6. 

185. 

6. 

190. 

7. 

195. 

4. 

200. 

3. 

205. 

3. 

210. 

8. 

215. 

1. 

220. 

6. 

225. 

0. 

230. 

5. 

235. 

1. 

240. 

7. 

245. 

1. 

260. 

2. 
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5.  MINIMIZING  A GENERAL  UNCONSTRAINED  FUNCTION  (OPTION  3) 

This  option  is  triggered  by  placing  UNCO  in  the  first  four  positions  of 
the  first  line  of  input.  The  main  difference  between  the  input  format  of  this 
and  that  for  an  indexed  function  (See  Section  4)  is  that  there  is  no  data 
defining  the  index  at  the  end  of  the  input. 

This  option  will  be  illustrated  after  a discussion  on  how  to  solve 
constrained  problems  by  solving  a sequence  of  unconstrained  problems. 


5.1  SOLVING  CONSTRAINED  PROBLEMS  VIA  UNCONSTRAINED  OPTIMIZATION 

Although  FACTNLS  is  primarily  for  unconstrained  problems,  it  is  possible 
to  use  it  as  an  aid  for  solving  constrained  optimization  problems.  There  is  a 
well-known  algorithm,  SUMT,  which  solves  constrained  problems  by  solving  a 
sequence  of  unconstrained  problems.  For  more  information  the  reader  is 
referred  to  [Fiacco  and  McCormick  1968]  or  [McCormick  1983,  Chapter  16].  A 
brief  description  of  the  SUMT  algorithm  follows. 

The  optimization  problem  is  assumed  to  be  in  the  form: 
minimize  f(x) 

^ s.t.  g^(x)  > 0,  for  i=l,  ...,  m, 

hj (x)  =0,  for  j^l , . . . , p. 

The  vector  x = , ...,  Any  optimization  problem  can  be  put  into  the 

above  form.  Denote  R°  = {x:  gi (x)  > 0,  for  i-1 , ...,  m} . The  SUMT  algorithm 
solves  the  problem 

m p 

minimize  P(x,rj^)  = f(x)  - rj^  E ln[g^(x)]  + E hj(x)/rj. 
x€R°  i=l  j=l 


for  a sequence  of  values  (rj^)  which  decrease  strictly  to  zero.  If  x(rj.) 
denotes  a solution  to  the  problem  for  rj^  , the  starting  point  for  the  new 
problem  with  r^^  + j^  is  x(rjj).  For  r^^  small,  x(rj.)  is  an  approximation  to  the 
solution  of  the  original  constrained  nonlinear  programming  problem.  There  is 
an  extrapolation  method  which  can  be  used  to  estimate  the  solution  accurately 
using  two  vectors  x(rjj)  and  x(rj.  + j^)  on  the  trajectory  of  unconstrained 
minimizers.  The  chemical  equilibrium  problem  will  be  used  to  explain  this. 

The  example  illustrating  this  option  is  to  solve  a constrained  problem 
by  solving  a sequence  of  unconstrained  problems.  The  theory  behind  this 
approach  is  given  after  a description  of  the  optimization  problem:  CHEMICAL 
EQUILIBRIUM.  The  example  and  data  are  taken  from  Bracken  and  McCormick 
[1968]  . 
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The  problem  of  determining  the  chemical  composition  of  a complex  mixture 
under  chemical  equilibrium  conditions  has  long  been  of  interest.  Such  problems 
arise  in  the  analysis  of  the  performance  of  fuels  and  propellants  and  in  the 
synthesis  of  complex  organic  compounds. 

A mixture  of  chemical  species  held  at  a constant  temperature  and 
pressure  reaches  its  chemical  equilibrium  state  concurrently  with  reduction  of 
the  free  energy  of  the  mixture  to  a minimum.  This  is  a consequence  of  the 
second  law  of  thermodynamics.  The  objective  function  to  be  minimized  in  the 
chemical  equilibrium  model  is  the  expression  of  the  free  energy  of  the 
chemical  mixture  under  study.  The  value  of  the  free  energy  of  the  mixture  is 
minimized  subject  to  the  chemical  reactions  possible  between  species  of  the 
mixture . 


White,  Johnson,  and  Dantzig  formulated  the  chemical  equilibrium  problem 
as  a mathematical  programming  problem  with  linear  mass  balance  constraints 
representing  the  possible  chemical  combinations  of  the  chemical  species  of  the 
mixture,  and  a nonlinear  objective  function  representing  the  free  energy  of 
the  mixture  (to  be  minimized) . They  investigated  steepest  descent  and 
piecewise  linear  programming  approaches  to  formulating  the  problem.  In  a 
second  paper,  they  explored  the  piecewise  linear  programming  problem  further. 
The  problem  is  discussed  briefly  by  Dantzig,  who  used  it  to  illustrate  the 
method  of  generalized  linear  programming. 


Consider  a mixture  of  m chemical  elements.  It  has  been  predetermined 
that  the  m different  types  of  atoms  can  combine  chemically  to  produce  n 
compounds,  where  the  monotonic  atom  is  regarded  for  our  purpose  as  a possible 
compound.  Define 


Xj  = the  number  of  moles  of  compound  j in  the  mixture  at  equilibrium, 
_ _ n 

X = the  total  number  of  moles  in  the  mixture,  where  x = S Xj  , 

j=l 

a^j  = the  number  of  atoms  of  element  i in  a molecule  of  compound  j , 
b^  = the  number  of  atomic  weights  of  element  i in  the  mixture. 


The  mass  balance  relationships  that  must  hold  among  the  m elements  are 


and 


n 

2 a^jXj  = bi  , i=l,  . . 

j=l 

. , m. 

(5.1) 

Xj  > 0,  j=l,  . 

. . , n. 

(5.2) 

Determination  of  the  composition  of  the  mixture  at  equilibrium  is 
equivalent  to  determination  of  the  values  of  Xj  (j=l,...,n)  that  satisfy  (5.1) 
and  (5.2)  and  also  minimize  the  total  free  energy  of  the  mixture.  The  total 
free  energy  of  the  mixture  is  given  by 


S Xj  [Cj  + in(Xj/x)]  (5.3) 

j=l 
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where 


Cj  = (F^/RT)  + in  P 

where  (Fj/RTj)  is  the  modal  standard  (Gibbs)  free  energy  function  for  the  j th 
compound,  which  may  be  found  in  tables,  and  P is  the  total  pressure  in 
atmospheres . 

Thus  the  nonlinear  programming  problem  is  as  follows.  Choose  Xj 
(j=l,...,n)  to  minimize  the  nonlinear  objective  function  (5.3)  subject  to 
linear  constraints  (5.1)  and  non-negativity  restrictions  (5.2). 

We  consider  the  example  problem  formulated  and  solved  by  White,  Johnson, 
and  Dantzig  [??].  We  solve  the  nonlinear  programming  problem  by  the 
sequential  unconstrained  minimization  technique  and  obtain  similar  answers. 

The  problem  considered  is  the  determination  of  the  equilibrium 
composition  resulting  from  subjecting  the  compound  I/2N2H4  + I/2O2  to  a 
temperature  of  3500“K  and  a pressure  of  750  psi.  In  Table  5 we  show  for  each 
compound  j of  10  possible  compounds  (where  the  monotonic  atoms  are  termed 
compounds)  the  Gibbs  free  energy  function  (F'’/RT)j  , the  computed  value  of  Cj 
for  P = 750  psi,  and  the  number  atoms  of  H,  N,  and  0 per  molecule.  The  number 
of  atomic  weights  of  H,  N,  and  0 in  the  mixture  are  bj^=2,  b2=l,  and  b3=l. 

Formulating  the  nonlinear  programming  model,  the  nonlinear  objective 
function  to  be  minimized  is 


Xj^  [ -6.089  + ln(xj^/x)] 


+ . . . 

+ Xio  [ -22.179  + ln(Xio/x)] 


and  the  linear  constraints  of  the  nonlinear  progranuning  problem  are  as 
follows : 


Xj^  + 2x2  + 2x3  + Xg  + Q = 2 , 


X4 

+ 2X5  + Xg  + Xy 

= 1. 

X3 

+ X7  + Xg  + 2X9  + Xio 

= 1. 

> 0 , X2  > 0 , . . . , x^  0 

> 0. 

Solving  the  above  nonlinear  programming  problem,  we  obtain  the  values  of 
Xj  (j=l,...,  10),  the  number  of  moles  of  the  10  compounds  present  in  the 
equilibrium  mixture,  which  are  given  in  Table  5.  These  values  agree  with 
those  obtained  in  [??].  The  corresponding  value  of  the  objective  function 
is  -47.76. 
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TABLE  5.  DATA  ON  1/2N2H^  + I/2O2  AT  3500"K,  750  psi. 


i=l 

i=2 

i=3 

j 

Compound 

(FVRT)j 

H 

N 

0 

1 

H 

-10.021 

- 6.089 

1 

0 

0 

2 

Hz 

-21.096 

-17.164 

2 

0 

0 

3 

H2O 

-37.986 

-34.054 

2 

0 

1 

4 

N 

- 9.846 

- 5.914 

0 

1 

0 

5 

Nz 

-28.653 

-24.721 

0 

2 

0 

6 

NH 

-18.918 

-14.986 

1 

1 

0 

7 

NO 

-28.032 

-24.100 

0 

1 

1 

8 

0 

-14.640 

-10.708 

0 

0 

1 

9 

O2 

-30.594 

-26.662 

0 

0 

2 

10 

OH 

-26.111 

-22.179 

1 

0 

1 

TABLE  6. 

COMPOSITION  OF  1/2N2H^  + I/2O2 
AT  3500’K,  750  PSI 


j 

Compound 

1 

H 

.0407 

2 

Hz 

.1477 

3 

H2O 

.7831 

4 

N 

.0014 

5 

Nz 

.4853 

6 

NH 

.0007 

7 

NO 

.0274 

8 

0 

.0180 

9 

O2 

.0373 

10 

OH 

.0969 

In  Figure  7 is  the  input  file  for  the  chemical  equilibrium  problem  with  an 
initial  value  of  rj^  = .01.  In  Appendix  A is  a portion  of  the  output.  The 
final  values  at  the  unconstrained  minimizer  are  then  used  in  the  input  file 
shown  in  Figure  8,  which  has  a value  of  "RK"  = .001.  The  solution  of  this 
second  unconstrained  problem  is  given  in  Figure  10.  With  these  two  vectors, 
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FIGURE  7.  INPUT  FILE  FOR  CHEMICAL  EQUILIBRIUM  PROBLEM  WITH  rk  = .01. 


UNCO --CHEMICAL  EQUILIBRIUM  PROBLEM  USING  BARRIER  FUNCTION 
SOURCE-  PP  253-254  OF  NONLINEAR  PROGRAMMING  BY  GPMCC 
C1*H+C2*H2+C3*H20+C4*N+C5*N2+C6*NH+C7*N0+C8*0+ 

C9*02+C 1 0*OH  +RK* ( BARLN (H) +BARLN ( H2 ) +BARLN ( H20 ) 

+BARLN(N)  + BARLN(N2)  + BARLN(NH)  + BARLN(NO)  + BARLN(O) 
+BARLN(02)  +BARLN(OH)  ) + ( (H+2 . *H2+2 . *H20+NH+0H-B1)**2 
+(N  +2.*N2  +NH  +NO-B2)**2  +(H20+N0+0+2 . *02+0H-B3)**2)/RK 
-ENTR0P(H+H2+H20+N+N2+NH+N0+0+02+0H)  + ENTROP(H)  + ENTR0P(H2) 

+ ENTR0P(H20)  + ENTROP(N)  + ENTROP(N2)  + ENTROP(NH)  + ENTROP(NO) 
+ ENTROP(O)  + ENTR0P(02)+  ENTROP(OH)  $ 

NUMBER  OF  UNKNOWNS,  CONSTANTS,  QUANTITIES,  AND  OBSERVATIONS 
10  14  0 0 

STARTING  VALUES  FOR  UNKNOWNS 


.1  . 

00001 

10. 

'H2' 

.1 

0.00001 

10. 

'H20' 

0.1 

0.00001 

10 

'N' 

.1 

0.00001 

10. 

'N2' 

0.1 

0.00001 

10 

'NH' 

.1 

0.00001 

10. 

'NO' 

.1 

0.00001 

10. 

'O' 

.1 

0.00001 

10. 

'02' 

0.1 

0.00001 

10 

'OH' 

.1 

0.00001 

10. 

CONSTANTS 


'Cl' 

-6.089 

'C2' 

-17.164 

'C3' 

-34.054 

'C4' 

-5.914 

'C5' 

-24.721 

'C6' 

-14.986 

'C7' 

-24.100 

'C8' 

-10.708 

'C9' 

-26.662 

'CIO' 

’ -22.179 

'Bl' 

2. 

'B2' 

1. 

'B3' 

1. 

'RK' 

.01 
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FIGURE  8:  INPUT  FILE  FOR  CHEMICAL  EQUILIBRIUM  PROBLEM  WITH  rk=.001 


UNCO- -CHEMICAL  EQUILIBRIUM  PROBLEM  USING  BARRIER  FUNCTION 
SOURCE-  PP  253-254  OF  NONLINEAR  PROGRAMMING  BY  GPMCC(RK=. 001) 
C1*H+C2*H2+C3*H20+C4*N+C5*N2+C6*NH+C7*N0=C8*0+ 

C9*02+C10*0H  +RK*(BARLN(H)+BARLN(H2)+BARLN(H20) 

+BARLN(N)  + BARLN(N2)+BARLN(NH)+BARLN(N0)+BARLN(0) 

+BARLN(02)  +BARLN(0H)  )+((H+2.  *H2+2 . *H20+NH+0H-B1)**2 
+(N  +2.  *N2  +NH  +N0-B2)**2  +(H20+N0+0+2 . *02+0h-B3)**2)/RK 
-ENTROP(H+H2+H20+N+N+N2+NH+N0+0+02+0H) 

+ ENTRO  P ( H ) +ENTROP ( N2 ) +ENTROP ( h2  0 ) +ENTROP ( N ) +ENTROP ( N2 ) +ENTROP ( NH ) +ENTROP 
ENTROP(O)  +ENTROP(02)+  ENTROP(OH)  $ 

NUMBER  OF  UNKNOWNS,  CONSTANTS,  QUANTITIES,  AND  OBSERVATIONS 
10  14  0 0 

STARTING  VALUES  FOR  UNKNOWNS 


'H' 

.049942  .00001  10. 

'H2' 

.15007  0.00001  10. 

'H20' 

.79116  0.00001  10. 

'N' 

.0066461  0.00001 

10 

'N2' 

.50722  0.00001 

10. 

'NH' 

.0050182  0.00001 

10 

'NO' 

.038836  0.00001 

10 

'0' 

.02846  0.00001 

10. 

'02' 

.052818  0.00001 

10 

'OH' 

.11177  0.00001 

10. 

CONSTANTS 

'Cl' 

-6.089 

'C2' 

-17.164 

'C3' 

-34.054 

'C4' 

-5.914 

'C5' 

-24.721 

'C6' 

-14.986 

'C7' 

-24.100 

'C8' 

-10.708 

'C9' 

-26.662 

'CIO' 

-22.179 

'Bl' 

2. 

'B2' 

1. 

'B3' 

1. 

'RK' 

.001 

VARIABLES 
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the  solution  can  be  approximated  using  extrapolation  theory.  The  two  point 
extrapolation  is 


x*(est)  = [cx(rj,  + i)  - x(ri^)]/(c  - 1) 

where 

c = + i . 

In  Table  7,  the  extrapolation  is  performed  using  the  output  from  the  two 
unconstrained  minimizations. 

Eventually,  solving  constrained  problems  will  be  an  automatic  procedure, 
For  now,  with  care,  the  FACTUNC  program  can  be  used  to  find  solutions  to 
constrained  problems. 

When  implementing  the  SUMT  approach,  it  is  important  to  use  the  single 
argument  function  BARLN(X) . This  prevents  the  sequence  of  iterates  from 
straying  from  R° , the  interior  of  the  feasible  region.  The  use  of  LN(X)  will 
result  in  disaster.  The  user  must  also  be  sure  that  the  initial  starting  point 
is  in  the  strict  interior  of  the  feasible  region. 

Choosing  the  and  + is  still  more  of  and  art  than  a science. 
Experience  is  often  a help  in  doing  this.  A large  value  to  start  with  is 
usually  better  then  a small  one,  which  creates  a number  of  difficult 
unconstrained  minimization  problems. 


6.  FUTURE  WORK 

The  system  described  in  this  manual  is  just  the  first  step  in  a proposed 
series  of  optimization  programs.  One  basic  addition  that  is  envisioned  is  the 
implementation  of  a general  indexing  capability  based  on  modern  data  base 
concepts.  Allowing  "subscripted”  data  will  extend  the  applicability  of  the 
system  to  large  problems  without  a concomitant  increase  in  program  size. 

The  second  generalization  is  to  the  solution  of  constrained  optimization 
problems.  Although  this  can  be  accomplished  (as  indicated  in  subsection  5.1) 
using  a barrier  function  approach,  it  is  important  for  a user  to  specify  in  a 
user  friendly  way  inequality  and  equality  constraints  on  the  problem  unknowns. 

A third  extension  will  be  to  allow  the  use  of  this  system  by  algorithmists  to 
test  their  methodology.  The  system  can  be  used  to  represent  the  problems  and 
to  compute  automatically  the  derivatives  required  by  the  algorithms. 
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TABLE  7. 

EXTRAPOLATION  FOR  CHENICAL  EQUILIBRIUM  PROBLEM 
USING  VALUES  FROM  TWO  SUMT  MINIMIZERS 


Compound 

Solution 
for  rij  = .01 

Solution 
for  rjj  = .001 

Extrapo- 

lation 

True 

Solution 

H 

.049942 

.041726 

.0408 

.0407 

Hz 

.15007 

.14829 

.1481 

.1477 

H2O 

.79116 

.78409 

.7833 

.7831 

N 

.0066461 

.0022263 

.0017 

.0014 

Nz 

.50722 

.48779 

.4856 

.4853 

NH 

.0050182 

.00008635 

-.0005 

.0007 

NO 

.038836 

.028593 

.0275 

.0274 

0 

.028460 

.019089 

.0180 

.0180 

O2 

.052818 

.038759 

.0372 

.0373 

OH 

.11177 

.098314 

.0968 

.0969 

Free 

Energy 

-48.702 

-47.850 

-47.755 

-47.76 
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APPENDIX  A: 


ABBREVIATED  OUTPUT  FOR  MAXL  OPTION  EXAMPLE. 


MAXL-FIND  PARAMETERS  DEFINING  FREQUENCY  FUNCTION  MIN  -(LOG  LIKE  FN) 
MINIMIZING  INDEXED  SUM 
HESSIAN  FORMED  EXPLICITLY 

FREQUENCE  FUNCTION  IS  WEIGHTED  SUM  OF  THREE  NORMALS 
(-FREQ)*LN((P/SIG)*NORDEN((Y-MU)/SIG)+(((l. -P-Q)/SIG2)* 

NORDEN( (Y-MU2)/SIG2) )+(Q/SIG3)*NORDEN( (Y-MU3)/SIG3) ) $ 

SINGLE  CHARACTER  OUTPUT  STRING  OF  (FUNCTION)  EQUATION 
( -FREQ)*LN( (P/SIG)*NORDEN( (Y-MU)/SIG)+( ( (1 . -P-Q)/S 
IG2)*NORDEN( (Y-MU2)/SIG2) )+(Q/SIG3)*NORDEN( (Y-MU3) 

/SIG3))$ 

NUMBER  OF  UNKNOWN  PARAMETERS , CONSTANTS , DATA  COLS.,  OBSERVATIONS 
THE  NUMBER  OF  UNKNOWNS  IS  = 8 


THE 

THE 

NUMBER  OF 
NUMBER  OF 

SYMBOLIC  CONSTANTS 
DATA  COLUMNS  = 2 

IN  THE  EQUATION 

(FUNCTION)  = 

THE 

NUMBER  OF 

DATA  OBSERVATIONS  = 

31 

PARAMETERS:  STARTING  VALUE,  LOWER 

AND  UPPER  BOUND 

MU 

.133243E+03 

.120000E+03 

.200000E+03 

SIG 

.947758E+01 

.800000E+01 

.700000E+02 

P 

.288020E+00 

.200000E+00 

.900000E+00 

MU2 

.171105E+03 

.141000E+03 

.230000E+03 

SIG2 

.337675E+02 

.lOOOOOE+02 

. lOOOOOE+03 

Q 

.200000E+00 

.lOOOOOE+00 

. 900000E+00 

MU3 

.200000E+03 

.lOOOOOE+03 

.300000E+03 

SIG3 

.300000E+02 

.lOOOOOE+02 

.500000E+02 

CONSTANTS  ( NONE  FOR  THIS  PROBLEM) 


INDEXED  DATA 
Y 

95.00000 

105.00000 

110.00000 

115.00000 

120.00000 

125.00000 

130.00000 

135.00000 

140.00000 

145.00000 

150.00000 

155.00000 

160.00000 

165.00000 

170.00000 


FREQ 

1.00000 

1.00000 

4.00000 

4.00000 

15.00000 

15.00000 

15.00000 

13.00000 

21.00000 
12.00000 

17.00000 

4.00000 

20.00000 

8.00000 
17.00000 
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175.00000 

8.00000 

180.00000 

6.00000 

185.00000 

6.00000 

190.00000 

7.00000 

195.00000 

4.00000 

200.00000 

3.00000 

205.00000 

3.00000 

210.00000 

8.00000 

215.00000 

1.00000 

220.00000 

6.00000 

225.00000 

.00000 

230.00000 

5.00000 

235.00000 

1.00000 

240.00000 

7.00000 

245.00000 

1.00000 

260.00000 

2.00000 

THE  OPTIMAL  VALUES  OF  THE  UNKNOWNS  ARE 

1 MU  = .12878E+03 

2 SIG  = .10297E+02 

3 P = .25475E+00 

4 MU2  = .15826E+03 

5 SIG2  = .21741E4-02 

6 Q = .15123E+00 

7 MU3  = .22224E+03 

8 SIG3  = .18505E+02 

THE  LEAST  SQUARES  VALUE  IS  = .11384E+04 

THE  GRADIENT  IS  -.37851E-13  -.60924E-13  -.60966E-11  -.43141E=13 

-.10113E-13 

THE  GRADIENT  IS  -.76739E-11  .67502E-13  -.20073E-12 
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APPENDIX  B: 


ABBREVIATED  OUTPUT  FOR  CHEKICAL  EQUILIBRIUM  PROBLEM  (rk=.01) 


UNCO- -CHEMICAL  EQUILIBRIUM  PROBLEM  USING  BARRIER  FUNCTION 
UNCONSTRAINED  PROBLEM 
HESSIAN  FORMED  EXPLICITLY 

SOURCE-  PP  253-254  OF  NONLINEAR  PROGRAMMING  BY  GPMCC 
Cl*H+C2*H2+C3*H20+C4*N+C5*N2+C6*NH-i-C7*N0+C8*0+ 

C9*O2+C10*OH  +RK*(BARLN(H)+BARLN(H2)+BARLN(H20) 

+BARLN(N)  + BARLN(N2)+BARLN(NH)+BARLN(N0)-»-BARLN(0) 

+BARLN(02)  +BARLN(OH)  )+( (H+2 .*H2+2 .*H20+NH+0H-B1)**2 
+(N  +2.*N2  +NH  +NO-B2)**2  +(H20+N0+0+2 .*02+0H-B3)**2)/RK 
- ENTROP ( H+H2+H20+N+N2+NH+N0+0-K)2+0H ) 

+ 

ENTROP ( H ) +ENTROP ( H2 ) +ENTROP ( H20 ) +ENTROP ( N ) +ENTROP ( N2 ) +ENTROP ( NH ) +ENTROP ( NO ) + 
ENTROP(O)  + ENTROP (02)+  ENTROP (OH)  $ 

SINGLE  CHARACTER  OUTPUT  STRING  OF  (FUNCTION)  EQUATION 
C1*H+C2*H2+C3*H20+C4*N+C5*N2+C6*NH+C7*N0+C8*0+C9*0 
2+C 1 0*OH+RK* ( BARLN ( H ) +BARLN ( H2 ) +BARLN ( H20 ) +BARLN ( N 
) +BARLN (N2 ) +BARLN ( NH ) +BARLN ( NO ) +BARLN ( 0 ) +BARLN ( 02 ) 

+BARLN(OH) )+( (H+2 . *H2+2 . *H20+NH+0H-B1 )**2+(N+2 . *N2 
+NH+N0-B2 ) **2+ (H20+N0+0+2 . *02+0H-B3 ) **2 ) /RK- ENTROP 
(H+H2+H20+N+N2+NH+N0+O+02+0H)+ENTR0P (H) +ENTROP (H2 ) 

+ENTROP ( H20 ) +ENTROP ( N ) +ENTROP ( N2 ) +ENTROP ( NH ) +ENTRO 
P ( NO ) +ENTROP ( 0 ) +ENTROP (02) +ENTROP ( OH ) $ 

NUMBER  OF  UNKNOWNS,  CONSTANTS,  QUANTITIES,  AND  OBSERVATIONS 


THE 

NUMBER  OF 

UNKNOWNS  IS  = 10 

THE 

NUMBER  OF 

SYMBOLIC  CONSTANTS  IN 

THE  EQUATION 

(FUNCTION)  = 

THE 

NUMBER  OF 

INDEPENDENT  VARIABLES 

= 0 

THE 

NUMBER  OF 

DATA  OBSERVATIONS  = 

0 

STARTING  VALUES  FOR  UNKNOWNS 

H 

.lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

H2 

.lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

H20 

. lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

N 

. lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

N2 

.lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

NH 

.lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

NO 

.lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

0 

.lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

02 

.lOOOOOE+00 

.lOOOOOE-04 

. lOOOOOE+02 

OH 

.lOOOOOE+00 

.lOOOOOE-04 

.lOOOOOE+02 

CONSTANTS 

Cl  -.608900E+01 
C2  -.171640E+02 
C3  -.340540E+02 
C4  -.591400E+01 


30 


C5  -.247210E+02 
C6  -.149860E+02 
C7  -.241000E+02 
C8  -.107080E+02 
C9  -.266620E+02 
CIO  -.221790E+02 
B1  .200000E+01 
B2  .lOOOOOE+01 
B3  .lOOOOOE+01 
RK  .lOOOOOE-01 


VARIABLES 

THE  OPTIMAL  VALUES  OF  THE  UNKNOWNS  ARE 


1 

H 

.49942E-01 

2 

H2 

.15007E+00 

3 

H20 

= 

.79116E+00 

4 

N 

.66461E-02 

5 

N2 

= 

.50722E+00 

6 

NH 

.50182E-02 

7 

NO 

— 

.38836E-01 

8 

0 

= 

.28460E-01 

9 

02 

.52818E-01 

10 

OH 

= 

.11177E+00 

THE  LEAST  SQUARES  VALUE  IS  = -.48702E+02 


THE  GRADIENT 
-.94860E-09 

IS 

-.22860E-07 

-.28291E-07 

-.22719E-08 

THE  GRADIENT 
-.49182E-07 

IS 

-.67860E-05 

-.49977E-09 

-.20685E-09 

-.10537E-03 

-.51097E-07 
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APPENDIX  C; 


PARTIAL  OUTPUT  FOR  CHEMICAL  EQUILIBRIUM  (rj^=  .001) 

UNCO- -CHEMICAL  EQUILIBRIUM  PROBLEM  USING  BARRIER  FUNCTION 
UNCONSTRAINED  PROBLEM 
HESSIAN  FORMED  EXPLICITYLY 

SOURCE-  PP  253-254  OF  NONLINEAR  PROGRAMMING  BY  GPMCC  (RK=.001) 
C1*H+C2*H2+C3*H20+C4*N+C5*N2+C6*NH+C7*N0+C8*0+ 

C9*O2+C10*OH  +RK*(BARLN(H)+BARLN(H2)+BARLN(H20) 

+BARLN(N)  + BARLN(N2)+BARLN(NH)+BARLN(N0)+BARLN(0) 

+BARLN(02)  +BARLN(OH)  )+((H+2.  *H2+2 . *H20+NH+0H-B1)**2 
+(N  +2.  *N2  +NH  +NO-B2)**2  +(H20+NCH-0+2 . *02+0H-B3)**2)/RK 
- ENTROP  (H+H2+H20+N+N2+NH+N0-K)+02-K)H) 

+ENTROP (H) +ENTROP (H2 ) +ENTROP (H20) +ENTROP (N) +ENTROP (N2 ) +ENTROP (N2 ) +ENTROP (NH) 
+ENTROP(NO)+ENTROP  (0)  + ENTROP  (02)+  ENTROP(OH)  $ 

SINGLE  CHARACTER  OUTPUT  STRING  OF  (FUNCTION)  EQUATION 
C1*H+C2*H2+C3*H20+C4*N+C5*N2+C6*NH+C7*N0+C8*0+C9*0 
2+C 1 0*OH+RK* ( BARLN ( H ) +BARLN ( H2 ) +BARLN ( H20 ) +BARLN (N 
) +BARLN ( N2 ) +BARLN ( NH ) +BARLN (NO ) +BARLN ( 0 ) +BARLN ( 02 ) 

+BARLN(OH) )+( (H+2 . *H2+2 . *H20+NH+0H-Bl)**2+(N+2 . *N2 
+NH+N0-B2)**2+(H20+N0+0+2 . *02+OH-B3)**2)/RK- ENTROP 
( H+H2+H20+N+N2+NH+N0+0+02+0H ) +ENTROP ( H ) +ENTROP ( H2 ) 

+ENTROP ( H20 ) +ENTROP ( N ) +ENTROP ( N2 ) +ENTROP (NH ) +ENTROP 
( NO ) +ENTROP ( 0 ) +ENTROP ( 02 ) +ENTROP ( OH ) $ 

NUMBER  OF  UNKNOWNS,  CONSTANTS,  QUANTITIES,  AND  OBSERVATIONS 
THE  NUMBER  OF  UNKNOWNS  IS  = 10 

THE  NUMBER  OF  SYMBOLIC  CONSTANTS  IN  THE  EQUATION  (FUNCTION)  = 14 
THE  NUMBER  OF  INDEPENDENT  VARIABLES  = 0 

THE  NUMBER  OF  DATA  OBSERVATIONS  = 0 

STARTING  VALUES  FOR  UNKNOWNS 


H 

.499420E-01 

.100000E-04 

.lOOOOOE+02 

H2 

.150070E+00 

.lOOOOOE-04 

.lOOOOOE+02 

H20 

.791160E+00 

.lOOOOOE-04 

.100000E+02 

N 

.664610E-02 

.lOOOOOE-04 

.100000E+02 

N2 

.507220E+00 

.lOOOOOE-04 

.lOOOOOE+02 

NH 

.501820E-02 

.lOOOOOE-04 

.100000E+02 

NO 

.388360E-01 

.100000E-04 

.100000E+02 

0 

.284600E-01 

.100000E-04 

.lOOOOOE+02 

02 

.528180E-01 

.100000E-04 

. 100000E+02 

OH 

.111770E+OO 

.100000E-04 

.100000E+02 

CONSTANTS 
Cl  -.608900E+01 
C2  -.171640E+02 
C3  -.340540E+02 
C4  -.591400E+01 
C5  -.247210E+02 
C6  -.149860E+02 
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C7  ~.241000E+02 
C8  -.107080E+02 
C9  =.266620E+02 
CIO  -.221790E+02 
B1  .200000E+01 
B2  .lOOOOOE+01 
B3  .lOOOOOE+01 
RK  .lOOOOOE-02 

THE  OPTIMAL  VALUES  OF  THE  UNKNOWNS  ARE 
THE  CURRENT  VALUES  OF  THE  UNKNOWS  ARE 


1 

H 

.41726E-01 

2 

H2 

.14829E+00 

3 

H20 

.78409E+00 

4 

N 

.22263E-02 

5 

N2 

= 

.48779E+00 

6 

NH 

= 

.86349E-04 

7 

NO 

= 

.28593E-01 

8 

0 

= 

.19089E-01 

9 

02 

.38759E-01 

10 

OH 

.98314E-01 

THE  LEAST  SQUARES  VALUE  IS  = -.47850E+02 


THE 

GRADIENT 

IS 

-.20583E-09 

-.21088E-08 

-.17735E-10 

-.20692 

THE 

GRADIENT 

IS 

-.13664E+02 

-.10834E-09 

-.24742E-08 

-.68331 
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